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As low-loss non-linear elements, Josephson junctions are the building blocks of superconducting 
qubits. The interaction of the qubit degree of freedom with the quasiparticles tunneling through 
the junction represent an intrinsic relaxation mechanism. We develop a general theory for the qubit 
decay rate induced by quasiparticles, and we study its dependence on the magnetic flux used to 
tune the qubit properties in devices such as the phase and flux qubits, the split transmon, and the 
fluxonium. Our estimates for the decay rate apply to both thermal equilibrium and non-equilibrium 
quasiparticles. We propose measuring the rate in a split transmon to obtain information on the 
possible non-equilibrium quasiparticle distribution. We also derive expressions for the shift in qubit 
frequency in the presence of quasiparticles. 
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I. INTRODUCTION 

The operability of a quantum device as a qubit re- 
quires long coherence times in comparison to the gate 
operation timeP^ Over the years, longer coherence times 
in superconducting qubits have been achieved by design- 
ing new systems in which the decoupling of the quantum 
oscillations of the order parameter from other low-energy 
degrees of freedom is enhanced. For example, in a trans- 
mon qubitP the sensitivity to background charge noise is 
suppressed relative to that of a Cooper pair box. Irre- 
spective of the particular design, in any superconduct- 
ing device the qubit degree of freedom can exchange en- 
ergy with quasiparticles. This intrinsic relaxation mech- 
anism is suppressed in thermal equilibrium at temper- 
atures much lower than the critical temperature, due 
to the exponential depletion of the quasiparticle popu- 
lation. However, both in superconducting qubits^ and 
resonators^ nonequilibrium quasiparticles have been ob- 
served which can lead to relaxation even at low tempera- 
tures. In this paper we study the quasiparticle relaxation 
mechanism in qubits based on Josephson junctions, both 
for equilibrium and nonequilibrium quasiparticles. 

Quasiparticle relaxation in a Cooper pair box was con- 
sidered in Ref. [51 In this system the charging energy is 
large compa red to the Josephson energy and quasiparti- 
cle poisoning^! is the elementary process of relaxation: 
a quasiparticle entering the Cooper pair box changes the 
parity (even or odd) of the state, bringing the qubit out of 
the computational space consisting of two charge states 
of the same parity. More recently the theory of Rcf. 5 
was extended to estimate the effect of quasiparticles in 
a transmonP In this case the dominant energy scale is 
the Josephson energy, so that quantum fluctuations of 
the phase are relatively small, while the uncertainty of 
charge in the qubit states is significant. As mentioned 
above, the advantage of the transmon is its low sensitiv- 
ity to charge noise. The possible role of nonequilibrium 
quasiparticles in superconducting qubits was investigated 
in Ref. 3; While the properties of many superconducting 
qubits - the phase and flux qubitsj^ the split transmon, 



and the newly developed fluxoniunP - can be tuned by 
an external magnetic flux, the effect of the latter on the 
quasiparticle relaxation rate has not been previously ana- 
lyzed. Elucidating the role of flux is the main goal of this 
work. In particular, we show that studying the flux de- 
pendence of the relaxation rate can provide information 
on the presence of nonequilibrium quasiparticles. 

The paper is organized as follows: in the next sec- 
tion we present results for the admittance of a Josephson 
junction and the general approach to calculate the de- 
cay rate and energy level shifts due to quasiparticles in 
a qubit with a single Josephson junction. In Sec. |III| we 
consider a weakly anharmonic qubit, such as phase qubit 
or transmon, and relate its decay rate, quality factor, and 
frequency shift to the admittance of the junction. The 
cases of a Cooper pair box (large charging energy) and of 
a flux qubit with large Josephson energy are examined in 
Sec. |IV| Some of the results presented in Sees. [n]|TV] have 
been reported previously^ in a brief format. In Sec. [v] 
we describe the generalization to multi-junction systems 
and study, as concrete examples, the two-junction split 
transmon and the many-junction fluxonium. We summa- 
rize the present work in Sec. |VI| Throughout the paper, 
we use units fi = ks = 1 (except otherwise noted). 



II. GENERAL THEORY FOR A 
SINGLE-JUNCTION QUBIT 

We consider a Josephson junction closed by an induc- 
tive loop, see Fig|T] The low -energy effective Hamilto- 
nian of the system can be separated into three parts 

H = H v + H q p + Ht ■ (1) 

The first term determines the dynamics of the phase de- 
gree of freedom in the absence of quasiparticles 

H V =4E C (n - -E. r cos£+^E L {0 - 2n<P e /<S> ) 2 , 

(2) 
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where N = —id/dip is the number operator of Cooper 
pairs passed across the junction, n g is the dimension- 
less gate voltage, $ e is the external flux threading the 
loop, $0 — h/2e is the flux quantum, and the parame- 
ters characterizing the qubit are the charging energy Eq, 
the Josephson energy Ej, and the inductive energy El- 
The second term in Eq. ([!]) is the sum of the BCS 
Hamiltonians for quasiparticles in the leads 



E 



qp ■ 



qp / j n ricr na ' 



(3) 
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=t, i accounts for spin. The quasiparti- 
(AJ) 2 , with^ and be- 



where d^^d^.) are quasiparticle annihilation (creation) 
operators and a 

cle energies are e J n 

ing the single-particle energy level n in the normal state 
of lead j, and the gap parameter in that lead, respec- 
tively. The occupations of the quasiparticle states are 
described by the distribution functions 



J=L,R, (4) 



assumed to be independent of spin; double angular 



brackets 



/qp 



denote averaging over the quasiparticle 



states. Hereinafter we assume for simplicity equal gaps 
in the leads, A L = A R = A. 

The last term in Eq. describes quasiparticle tunnel- 
ing across the junction and couples the phase and quasi- 
particle degrees of freedom 
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(5) 



The electron tunneling amplitude t in this equation de- 
termines the junction conductance, gx — Aire 2 v L v R t 2 in 
the tunneling limit f< 1 which we are considering. From 
now on, we assume identical densities of states per spin 
direction in the leads, v L — v r — vq. The Bogoliubov am- 
plitudes v? rL , v 3 n can be taken real, since Eq. ^ already 
accounts explicitly for the phases of the order parameters 
in the leads via the gauge-invariant phase difference^ in 
the exponentials. Accounting for the Josephson effect 
and quasiparticles dynamics by Eqs. ([2|-([5]) is possible 
as long as the qubit energy ui and characteristic energy 
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FIG. 1: (a) Schematic representation of a qubit controlled by 
a magnetic flux, see Eq. pi), (b) Effective circuit diagram 
with three parallel elements - capacitor, Josephson junction, 
and inductor - characterized by their respective admittances. 



5E of quasiparticles (as determined by their distribution 
function and measured from A) are small compared to 
AP oj, 6E <C 2A. In this low-energy limit, we may fur- 
ther approximate u 3 m ~ v 3 n ~ \j\j2. Then the operators 
e ±%( pl 2 in Eq. ([H]), which describe transfer of charge ±e 
across the junction, combine to give 



na ma 



H.c. 



(6) 



Starting from this low-energy tunneling Hamiltonian, in 
the next section we calculate the dissipative part of the 
junction admittance. 



A. Response to a classical time-dependent phase 

We consider here the "classical" dissipative response 
of a Josephson junction to a small ac bias to show that 
Eq. |6j correctly accounts for the knowrP^junction losses 
in the low-energy regime. These "classical" losses are 
directly related to the decay rate in the quantum regime, 
as we explicitly show in the next section. 

We assume a time-dependent bias v(t) = vcos(ujt) of 
frequency uj > superimposed to a fixed phase differ- 
ence i^o ■ in other words, we take the phase to be a 
time-dependent number which, by the Josephson equa- 
tion dip/dt — 2ev(t), has the form 



ip{t) = (po H sin(ti)t) 



(7) 



Here we focus on the linear in v response in the low- 
energy regime. Expressions for the current through the 
junction valid beyond linear response can be found, for 
example, in Ref. [TT] Substituting Eq. ([7| into Eq. (RjJ), 
expanding for small v , and keeping the linear term, we 
find for the time-dependent perturbation 6H(t) causing 
the dissipation 



SH(t) 

Hac = a cos 



H AC sin(wi) , 
2 oj 



E 



a L ^a R +Hc 



(8) 



The average dissipated power can be calculated using 
Fermi's golden rule: it is given by the product of the 
transition rate times the energy change in a transition 
between quasiparticle states caused by the perturbation. 
The energy change in a transition is ±w by energy con- 
servation, with the two signs corresponds to the events 
giving energy to or taking energy from the system. The 
average power P is 



P = 2nJ2((\(W qP \H A c\{v} 



qp/ 



OJ 



(9) 



A,qp 



E n 



where E v qp and £\,qp are the total energies of the quasi- 
particles in their respective initial {?y} qp and final {A} qp 
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states. We use Eq. |8| to evaluate the matrix element, 
average over initial quasiparticle states, and sum over fi- 
nal states to find 



P 



1 



ReYj(u],<p )v 2 



witrP 



1 + cos tp 



ReF qp (u;) 



(10) 



(11) 



Here Re Y qp is the real part of the quasiparticle contribu- 
tion to the junction admittance at zero phase difference, 



Re Y qp (w) =g T 



2A 



dx 



1 



o yiyi + w/A 

[f E ((l + x)A)-f E ((l + x + oj/A)A)]. 

(12) 

In deriving these formulas we have approximated the 
standard BCS density of states functions as 



A 2 ' V? 




1 

2x 



and taken equal quasiparticle occupations in the two 
leads, f L = f = /; we use this simplifying assumption 
throughout the paper. We indicate with f E the energy 
mode of the distribution function 



Me) = i[/(0 + /(-£)], 



(14) 



where e = ^/£ 2 + A 2 . Equation jllj for the real part of 
the admittance, valid at lo > 0, agrees with the linear 
response, low-energy limit of the non-linear I-V charac- 
teristic presented in Ref. ITj Extension to u < is found 
by noticing that Re Y qp is an even function of frequency. 

In thermal equilibrium and at low temperatures T <C 
A the distribution function can be approximated as 



f E (e) 



and Eq. (12) gives, at arbitrary ratio lj/T 
ReY^)=g T 2 ^e-^e^K Q (M, 



(15) 

■j/T _ 

(16) 



Here Kq is the modified Bessel function of the second 
kind with asymptotes 



K (x) 



e x v / tt/2x , x > 1 
ln2/x — 7_e , i<1 



with 7s the Euler gamma. 

For a generic distribution function, we can relate 
Re Yqp to the density of quasiparticle n qp in the high- 
frequency regime u 3> SE, where SE indicates the char- 
acteristic energy of quasiparticle (measured from the gap) 



above which the occupation of the quasiparticle states 
can be neglected; in thermal equilibrium SE ~ T. Under 



the assumption u) 3> SE we obtain from Eq. (12 1 



Rey#(w) * 



2A 



where 



•-qp 



2v A 



3/2 



(18) 



(19) 



is the quasiparticle density normalized to the Cooper pair 
density and 



2V2i/ A 



dx 



f E ((l + x)A) 



(20) 



is the density written using the approximation in 
Eq. (13). Note that in thermal equilibrium at low tem- 



peratures, Eq. (15), we have 



n qp = 2v V2irATe- A/T . 



(21) 



Then using Eq. (17), it is easy to check that for T <C u> 



(13) Eq. (fl6|) takes the form given in Eq. (fl8|) 



The real and imaginary parts of the admittance satisfy 
the Kramers-Kronig relations. However, when taking the 
Kramcrs-Kronig transform of the real part, a purely in- 
ductive contribution to the imaginary part can be missed. 
Indeed, at low energies the complex junction admittance 
(obtained from the expressions in Ref. Ill p can be written 



1 



1P COS f + Yqp(w) — 



- cos If 



VjJLj 



where 



' qp 



= !e{A) 



(22) 



(23) 



can be interpreted as the population of the Andreev 
bound states^ and the inverse of the Josephson induc- 
tance is 

J- = . 9 T^A qp (24) 

(the subscript qp in A qp is used to indicate that in this 
expression it may be necessary to account for the effect 
of quasiparticles on the gap, see Sees. II C and IIIB). 



Unlike the Andreev states, free quasiparticles con- 
tribute to both dissipative and non-dissipative parts 
of the total admittance Yj via the complex term Y qp . 
The real part of the quasiparticle admittance is defined 
in Eq. (12), while its imaginary part is given by the 



(17) Kramers-Kronig transform of that expression, 



ImFqp(w) 



2AP 

-9t 

-f E ((l + y)A) 



f E ((l + x)A) 



dx f 00 dy 

i i 

x — y + lu/A x — y 



(25) 
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where P denotes the principal part and uj > 0. Using that 
Im Yq P is an odd function of frequency, we can simplify 
the above expression to a form with a single rather than 
double integral 



lmY qp (uj) = g T - 



2A 



M/ * f E ((l + x)A) 



dx 



:^\ujJ/A 



qp 



(26) 

As discussed above for the real part, an analytic ex- 
pression for ImF q p can be obtained in thermal equilib- 
rium, 



ImY q e p >) = ~9T 



2A 



UJ 



-A/T 

e 'it 



1 



-M/2T /o 



H 
2T 



(27) 

Here Jo is the modified Bessel function of the first kind 
with asymptotes 



( , J e x y/l/2irx, x> 1 
[X, -\l + x 2 /A, x<l 



(28) 



For arbitrary distribution function satisfying the high- 
frequency condition uj 3> SE we find 



hnY»(w 



1 



2A 



~9t — 
2 uj 




2^qp 



(29) 



Using Eq. (21 1 and the large- a; limit in Eq. (28), it is easy 



to show that for T <C to Eq. (27) reduces to the general 
expression in Eq. ( 29 ). In the high- frequency regime, real 



and imaginary parts of the quasiparticle admittance can 
be combined into the complex admittance 



iujL , 



(30) 



By substituting Eq. (I30J) into Eq. (22) we find that in 



the total admittance Yj the coefficient multiplying x qp is 
proportional to (1 — cost/?) and vanishes for tp = 0. This is 
in agreement with the absence of Andreev bound states 
when there is no phase difference across the junction. 



B. Transition rates 

The effects of the interaction between quasiparticles 
and qubit degree of freedom, Eq. ([5| , can be treated per- 
turbativcly in the tunneling amplitude t. The interaction 
makes possible, for example, a transition between two 
qubit states (initial, and final, |/), differing in energy 
by amount ujif > 0) by exciting a quasiparticle during 
a tunneling event. The rate for the transition between 
qubit states can be calculated using Fermi's golden rule 



fqp/ 



r^ /= 27T J2 ((|(/,{A}qp|ffT|*,Mc 

{A} qP 

X § (E\,qp — ^?7,qp — W l/)))qp • 



(31) 



We remind that in our notation E Vtqp (E\ iqp ) is the total 
energy of the quasiparticles in their initial (final) state 
{^Iqp ({^}qp); an d double angular brackets ((. . .)) qp de- 
note averaging over the initial quasiparticle states whose 
occupation is determined by the distribution function. 

In the low-energy regime we are considering, the tran- 
sition rate factorizes into terms accounting separately for 
qubit dynamic and quasiparticle kinetics 



(/|sm||i) 



S qp (w if ) 



(32) 



Equation ( 32 ) is one of the main results of this work: it 



shows that the qubit properties affect the transition rate 
via the wavefunctions \i), and |/) entering the matrix 
element, while the quasiparticle kinetics is accounted for 
by the quasiparticle current spectral density S'qp 



Sqp M 



16Ej 



1 



oo 

dx — 

II y/x^X + UJ / A 



f E ((l + x)A) 



x (l-/ B ((l + aOA+w)) 



where uj > and we used the relation 



Ej = g T A/8g K 



(33) 



(34) 



with = e 2 /2ir the conductance quantum. The expres- 
sion for S'qp at uj < is obtained by the replacements 
x — > x — uj/A, uj — s- —uj in the integrand in Eq. (33 ). 

The spectral density S qp depends on the detail of the 
distribution functions. In thermal equilibrium at low 



temperatures T <C A, using Eq. (15) we find 



Note that the equality 

SqpH 



(35) 



(36) 



implies that in thermal equilibrium the transition rates 
are related by detailed balance, 



/- 



(37) 



The similarity between Eq. (35) for S'qp and Eq. (16) 
for ReYqp is not accidental. In thermal equilibrium the 



following fluctuation-dissipation relation holds 
uj 1 

7T 9K 



S e q l(u) + S e q l(-uj) = ^—ReY^(uj)coth (^) . (38) 



Moreover, in the low-energy regime for an arbitrary dis- 
tribution function the two quantities are also related by 



uj 1 

1" 9K 



S qp (uj) - S qp (-uj) = Rer qP M . (39) 



5 



In the high-frequency regime ui ^> SE, we can simplify 
the above relation to 



4P 7T iP 



(40) 



For the transition rates this corresponds to neglecting the 
downward transitions with w,/ < 0, in which a quasipar- 
ticle looses energy to the qubit, compared to the upward 
ones. This is a good approximation since the assumption 
uj 3> SE means that there are no quasiparticles with en- 



ergy high enough to excite the qubit. Equation (40) can 
be checked by comparing Eq. (181 to 



(41) 



with Ej given in Eq. ( 34 ) and the the normalized quasi- 
particle density x qp in Eq. (19). 



C. Energy level corrections 

In addition to causing transitions between qubit lev- 
els, the quasiparticles affect the energy E^ of each level 
i of the system. We can distinguish two quasiparticle 
mechanisms that modify the qubit spectrum and hence 
separate two terms in the correction SEi to the energy, 



SEi = 5E i>Ej + SE, 



«,qp • 



(42) 



First, in the presence of quasiparticles the Josephson 
energy takes the form 



E 



J,qp 



9T 



8fJ 



K 



A qP (l 



2xt 



(43) 



with defined in Eq. (|23|). As mentioned after Eq. (j24h , 



we use A qp to distinguish the self-consistent gap in the 
presence of quasiparticles from the gap A when there are 
no quasiparticles. At leading order in the quasiparticle 
density we have 



A(l-Zqp) 



(44) 



Treating these modifications to the Josephson energy as 
perturbations, the correction to the energy of level i is 



SE i}Ej = Ej (a 



2xi 



{l\ cos (p\i) 



(45) 



Second, the virtual transitions between the qubit levels 
mediated by quasiparticle tunneling cause a correction 
that can be expressed in terms of the matrix elements of 
sin 0/2 as 



SE, 



«,qp 



E 

k^i 



where 



Wife 



I sin 



Ek — Ei 



(46) 



(47) 



The derivation of the above formulas and the definition 
of function F qp in terms of the quasiparticle distribution 



function [Eq. (A19)] are given in Appendix [Aj Here we 
give the relation between F qp and the imaginary part of 
the quasiparticle impedance, 

F qp (u)+F qp (-Lu) = -^—ImY qp (oj), (48) 

which we will use in the next section to obtain the 
quasiparticle-induced change in the qubit frequency. 



III. SINGLE JUNCTION: WEAKLY 
ANHARMONIC QUBIT 

As an application of the general approach described 
in the previous section, we consider here a weakly an- 
harmonic qubit, such as the transmon and phase qubits. 
We start with the the semiclassical limit, i.e., we assume 
that the potential energy terms in Eq. ([2| dominate the 
kinetic energy term proportional to Eq- This limit al- 
ready reveals a non-trivial dependence of relaxation on 
flux. Note that assuming El i= we can eliminate n g 
in Eq. |2j) by a gauge transformation.^ I n the transmon 
we have El — and the spectrum depends on n g , dis- 
playing both well separated and nearly degenerate states, 
see Fig. [2] The results of this section can be applied to 
the single-junction transmon when considering well sep- 
arated states. The transition rate between these states 
and the corresponding frequency shift are dependent on 
rig. However, since Ec <C Ej this dependence introduces 
only small corrections to r„_>„_i and Slj; the corrections 
are exponential in —y/8Ej/Ec- By contrast, the lead- 
ing term in the rate of transitions T e ^ a between the even 
and odd states is exponentially small. The rate r e .<_>. of 
parity switching is discussed in detail in Appendix [C| 

The potential energy in Eq. ^ is extremized at phase 
ipo satisfying 



Ej sin ipa + E L (ip - 27r$ e /$ ) = 0. 



(49) 



For Ej < El there is only one solution at the global 
minimum. For Ej > El however, there can be mul- 
tiple minima; their number depends both on the ratio 
Ej/El and the external flux <I> e . Here we assume that 
the flux is such that distinct minima are not degenerate; 
in particular, this means that the flux is tuned away from 
odd integer multiples of half the flux quantum.^ For the 
transmon with El = 0, we can take = as solution to 
Eq. (49). Next, we expand the potential energy around 



a minimum and find at quadratic order 

= AE c n 2 + \{El + Ej cos <po) (0 ~ <p ) 2 . (50) 

Fluctuations of the phase around ipo are small under the 
assumption 



E c 
n— < 1, 

Wio 



(51) 
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FIG. 2: Schematic representation of the transmon low en- 
ergy spectrum as function of the dimensionless gate voltage 
n g . Solid (dashed) lines denotes even (odd) states (see also 



Sec. IV A I. The amplitudes of the oscillations of the energy 
levels are exponentially smallP see Appendix [b] here they 
are enhanced for clarity. Quasiparticle tunneling changes the 
parity of the qubit sate. The results of Sec. |III| are valid for 
transitions between states separated by energy of the order of 
the plasma frequency lo p , Eq. ( |56[ ), and give, for example, the 
rate Fi_>o. For the transition rates between nearly degenerate 
states of opposite parity, such cLS To — Ye.} S66 Appendix [(5| 



where n denotes the energy level and 

u w = y/8E c {El + Ej cosp ) 



(52) 



is the qubit frequency in the harmonic approximation. 
Note that anharmonicity and quality factor Q determine 
the operability of the system as a qubitPThe anharmonic 
correction to the transition frequencies can be calculated 
by considering the effect on the spectrum of the next or- 
der in the expansion around po (cubic for the phase qubit, 
quartic for the transmon), which defines an anharmonic 
potential well of finite depth U. Then the operability 
condition can be expressed as Q/n w 3> 1, where n w is 
the number of states in the potential well, n w ~ ?7/u;xo- 
In a weakly anharmonic system, n w can be large; how- 
ever, if the quality factor is larger the system can be used 
as a qubit despite the weak anharmonicity, as it is indeed 
the case for the transmonP 



The condition for small phase fluctuations in Eq. (51) 
enables us to calculate the matrix element of operator 
smp/2 by expanding around tpo up to the second order 
and using standard expressions for the matrix elements 
of the position operator between eigenstates |n), |m) of 
the harmonic oscillator [cf. Eq. ( |50|]. To first order in 
Ec/(j->w we find (see also Appendix |D[) 



(m| sin — \n) 



2 






u w \ 2) 



1 



COS (fo 



, Ec r x r / LlU ,1 + cos <A) 
H [ndm^n-i + (n + 1) o m , n +iJ ~ ■ 

U>1Q & 

(53) 

Note that in the first term on the right hand side the 



corrections due to the non-linearity of sine (the second 
term inside the square brackets) are indeed small if condi- 
tion (51 1 is satisfied. In addition, we have neglected here 
the anharmonic corrections to the states used to calcu- 
late the matrix element; this is a good approximation for 
low-lying levels n -C n w W^ For the transmon (po — 0) the 
leading term in Eq. ( 53 ) is of linear order in Ec /wio ^ lj 



as we show in Appendix [E] by including the first anhar- 
monic correction to the states, the next non-vanishing 
term in the square of the transmon matrix element is cu- 
bic in Ec/uxo, rather than quadratic as for the harmonic 
oscillator. Therefore in the case of the transmon keep- 
ing only the leading term is a better approximation than 
naively expected. 

Equation ( 53 ) shows that at leading order we can re- 



strict our attention to transitions involving only neigh- 
boring levels. Concentrating here on low-lying levels, us- 
ing Eqs. ( 32 ) , ( 39 1 , and ( 53 1 we find the following relation 
between transition rate and impedance 



n 1 
-ReF qp (wi )- 



cos tpo 



(54) 



where we also used Ec = e 2 /2C. In the high-frequency 
regime, the upward transition rate can be neglected, 
r n _x_> n — 0, and the above expression simplifies to^ 
[see also Eq. ( 40 ) and the text that follows it] 



r, 



^Re^/Ko)- 



cos ipo 




(55) 



WlO 2-7T 



(1 + COSt^o) 



In the last expression we used Eq. (18) and introduced 
the plasma frequency 



\/8EcEj 



(56) 



The above equation can also be obtained by substituting 
directly Eq. ((401) into Eq. p2l). For n = 1 and p = 



Eq. ( 55 ) reduces to the transition rate presented in Ref. [31 



The transition rate in Eq. ( 55 1 is proportional to the 
(possibly non-equilibrium) quasiparticle density a; qp and 
depends on the external flux $ e via (po and uiiot see 
Eqs. (49) and (52). The flux dependence is in general 
sensitive to the states involved in the transition. This 
sensitivity can already be seen for transitions between 
harmonic oscillator states: due to the non-linear interac- 
tion between phase and quasiparticles, see Eq. ([6]), tran- 
sitions between distant levels are possible. These tran- 
sitions are suppressed by the smallness of phase fluctua- 
tions when Ec/ujiq <C 1. For example, the rate for the 
2^0 transition is 



r 2 



7T 9k 4P Wo 



1 



COS po 



(57) 



Note that in contrast to Eq. (55), Eq. (57) cannot be 
written in terms of the real part of the total admittance 
of the junction: while in Eq. (55) the phase enters via 



7 



the factor (1 -fcos^o) as in Eq. (Ill, in Eq. (57) ReF qp 
is multiplied by (1 — cost^o)- To obtain T2_>.o we sub- 
stituted into Eq. (|32l the high-frequency relation (|40|, 



while the explicit form of the squared matrix element 
|(0| sin(<£/2)|2)| 2 is found by setting n = 2, and keeping 
the leading term in Ec/^io, in the formula 





2 


(Ec\ 


(0|sin||n) 


— e "io 











1 — ( — 1)™ costpo 
2n\ 



(58) 



derived in Appendix [Dj Equation ( 58 ) is valid for any 
ratio Ec/uiq for transitions between eigenstates of the 
harmonic oscillator. When (fo = 0, Eq. (58) gives van- 
ishing matrix elements for even n - this is an example of 
the more general selection rule according to which only 
transitions between states of different parity are allowed 
at ipo = 0. 

The rate for transitions between excited states and the 
ground state in the case of large phase fluctuations can 
be obtained using Eq. (58) when E L <C E c and Ej < 
lulc — \/&EcEl- The latter condition enables us to 
neglect the Josephson energy term in Eq. (|2j). Then using 
Eqs. ( |32[ ) and ( |41[ ) with cjjj = nuj^c we find that the 
transition rate has a maximum for n = n with n ~ 




-(-l)"cos27r$ e /$ ] 
(n-n ) 2 ' 



2n n 



(59) 



Here we have approximated e y y n /n\ 



exp[— (n 



y) /2j/]/v27r y; the appr oxim ation is valid for y 3> 1 and 



-y 



< 



2y. Equation (59) shows that when the charg- 



ing energy is the dominant energy scale, dissipation is the 
strongest for transitions between states whose energy dif- 
ference (nuiLc) corresponds to the energy change (Ec) 
caused by the transfer of a single electron through the 
barrier, as in the "quasiparticle poisoning" picture for 
the Cooper pair boxP We stress that in the present case 
charge is not quantized, due to the finite value of the in- 
ductive energy -E^P^ We will comment on the relation 
between Eq. (59) and the transition rate in the Cooper 



pair box in Sec. IV A 



A. Quality factor 

Returning now to the semiclassical regime of small Eq, 
Eq. ( 55 ) with n = 1 enables us to evaluate, in the high- 



frequency regime, the inverse Q-factor for the transition 
between the qubit states 



1 

Qio 



rv 



l 



ReY^( 



E r: 1 



COS Ifo 



(60) 



We stress that this formula is valid not only in thermal 
equilibrium, but also in the presence of non-equilibrium 



quasiparticles with characteristic energy 8E <C wiq. We 



can generalize Eq. ( 60 1 to account for the possible coexis- 
tence of non-equilibrium and thermal quasiparticles. We 
take the distribution function in the form 



f E (e) = f ne (e) + f eq (e) . 



(61) 



where f ne is the non-equilibrium contribution, insensitive 
to temperature and satisfying the high-frequency condi- 
tion ujiq 3> 5E, and f eq is the equilibrium distribution 
of Eq. (15). Noting that within our assumption the two 



terms in Je contribute separately to the transition rates 
and that for the thermal part we cannot in general neglect 
the "upward" transitions, using Eqs. (32), (35), (41), and 



(53) we find 



1 



r - 



1 + COS Ifo w p 



^10 



2tt 



-10 




(62) 



where x ne is the normalized non-equilibrium quasiparti- 
cle density [cf. Eq. (T9J)]. 

Recently good agreement between theory, Eq 



(62), 



and experiment has been shown for single-junction trans- 
mons (ipo = 0, wio = oj p ) in the temperature range 
10-210 mK.SSI However, while these measurements indi- 
cate that thermal quasiparticles are the main cause of 
relaxation above ~ 150 mK, one cannot conclude that 
non-equilibrium quasiparticles are present from the lower 
temperature data: by Matthiessen rule, any other re- 
laxation mechanism which is independent of (or weakly 
dependent on) temperature would have the same limit- 
ing effect on Qio as the first term in square brackets in 
Eq, " 



(62) 



VA 



As we will discuss in more detail in Sec 
similar measurements on a flux-sensitive device shouk 
enable one to decide on the presence of non-equilibrium 
quasiparticles, since Eq. (62) [and its analogous for the 
split transmon, Eq. ( 127)]describes the effect of flux on 



both equilibrium and non-equilibrium quasiparticle con- 
tributions to Qio, and other sources of relaxation respond 
differently to the flux. 



B. Frequency shift 

A further test of the theory presented in Sec. [TT] is 
provided by the measurement of the qubit resonant fre- 
quency. In the semiclassical regime of small Ec, the qubit 
can be described by the effective circuit of Fig. [ijb) , with 
the junction admittance Yj of Eq. (22), Yc = iojC, and 
Yi = 1/iojL [the inductance is related to the inductive 
energy by E L = ($ /27r) 2 /L]. As discussed in Ref. [TO1 



for parallel elements the total admittance Y is the sum 
of their admittances, 



Y = Yj + Y c + Y L 



(63) 



and the resonant frequency uj r is the zero of the total 
admittance, Y(u r ) = 0. In the absence of quasiparticles 
we find uj r = u>\a with wio of Eq. (52). 



In the presence of quasiparticles, by considering their 
effect on the junction admittance at linear order in the 
quasiparticle density £ qp and Andreev level occupation 



x qp we obtain 



Soj 



with 



Y 

2C qp 



(wio) 



I + cos ipo 7rgrA 
2 Clo w ' 

_ 7rgrA 
2Cw 10 



, cos (fio 



x qp COS if 



(64) 



(65) 



The last term in Eq. ( 65 ) originates from the gap suppres- 
sion by quasiparticlesfcf. Eq. (44)]. This term was ne- 
glected in Ref.[in]as it is subleading in the high-frequency 
regime considered there [see Eq. ( 73 )] . The correction Slo 
has both real and imaginary parts. The imagina ry p art 
coincides^ with half the dissipation rate in Eq. ( 55 1 for 
the n = 1 — > transition. Here we show that the real part 
of 8uj r obtained in the effective circuit approach agrees 
with the quantum mechanical calculation. 

Within the harmonic approximation of Eq. ( 50 ) , the 



energy difference uit between the neighboring levels -Ej+i 
and Ea. 



For the quasiparticle tunneling term, we substitute 
Eq. ([53} into Eq. ((46} to get 



Su)i 



Ec 



[F q p(wio) + Fqp(-UJlo)] 



1 



COS (fo 



(71) 



Finally, using the relation ( 48 ) and adding the two terms 
we arrive at 



x 1 t v ( Ni + cos^o 

= ~ 2c 1 qp( Wl °) 2 

- -z — COS (p [x 
2 wio 



(72) 



MP 1 ^^qp) 



This expression agrees with the real part of Eq. ( 65 ) . We 



note that by extending the above consideration to include 
the next order in Ec /ojio, anharmonic corrections to the 
spectrum can be calculated. They are dominated by the 
anharmonicity of the cosine potential in Eq. ([2]), with 
quasiparticles contributing negligible additional correc- 
tions. For the case of the transmon, the leading anhar- 
monicity can be found in Ref. O 



In the high-frequency regime, using Eq. ( 29 ) the rela- 
tive frequency shift is 



= Ei + i — Ei 



(66) 



is of course independent of the level index i. The quasi- 
particle corrections to energy levels of Sec. |II C| cause a 
correction Soji to uii, 



Sljj = 5E, 



i+i 



5E X . 



(67) 



As we show below, at leading order in Eq/ojiq this cor- 
rection is also independent of level index, i.e, it represents 
a renormalization of the system resonant frequency. 
As in Eq. (42 1, we separate the contributions due to 



change in the Josephson energy and due to quasiparticle 
tunneling, 



u *,qp- 



(68) 



For the first term on the right hand side, we use Eq. (45) 



together with the matrix element of cos <p at first order 
in E c /u w [see Eq. flD9}], 



(i\ cos <p\i) ~ cos <pa 



1 



to find 



5ui 



1 2 
1 



-Ec 



apo (a 



2a; qp) 



(69) 



(70) 



As discussed in Sec. II C the term proportional to x qp is 
due to the ga_ 



d suppression in the presence of quasipar- 
ticles, Eq. (|44), while x^ p accounts for the occupation of 



the AndreevBound states. 



15 
'24 



1 qp 



(1 - COS^q) 



1 + COS ipo 

2tt 



COS ifio 



(73) 



Note that in the limit ojiq < Awe can neglect the cosine 
compared to the term multiplied by square root inside 
round brackets. However, this cosine term is the appro- 
priate subleading contribution, since the terms neglected 
in deriving the energy corrections presented in Sec. |II C| 
are suppressed by luiq / A with respect to the leading con- 
tribution. 

In recent experiments with single-junction trans- 
mon^22l relative shifts of order 10 -5 have been measured 
at temperatures ~ 200 mK, in agreement with Eq. (72). 



Together with the above mentioned measurements of the 
transition rates in the same devices, this is an additional, 
independent check of the validity of the present theory in 
the regime T > 150 mK. While in the transmon (ipo = 0) 
there are no Andreev bound states [indeed, in this case 
their contribution to the frequency shift is absent, see 
Eq. (73)], in a phase qubit both Andreev levels occu- 
pation Xq p and free quasiparticle density x qp affect the 
frequency. Assuming that the two quantity are propor- 



tional, 



,the ratio between frequency shift, 



Eq. (73), and transition rate, Eq. (55), in the high fre- 



quency regime is independent of the quasiparticle density. 
The constancy of this ratio has been recently verified by 
injecting a variable (but unknown) number quasiparticles 
in a phase qubit piJ 
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IV. SINGLE JUNCTION: STRONG 
ANHARMONICITY 

Here we consider the regime, complementary to that 
of the previous section, of qubits with large anharmonic- 
ities. We study first the single junction Cooper pair box 
(CPB); as for the transmon, it is insensitive to flux, but in 
contrast to the transmon the CPB properties are strongly 
affected by the value of the dimensionless gate voltage n g . 
Then we analyze a flux qubit, for which the external flux 
is tuned near half the flux quantum, $ e s» $o/2- 
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A. Cooper pair box 

The CPB is described by Eq. ^ with E L = and 
Ec ^> Ej. In this limit, it is convenient to rewrite the 
Hamiltonian in the charge basis as^ 



H=E c Y / (q-2n g ) 2 \q)(q\ 
q 

-^£(k)<? + 2| + k + 2><<z|). 



(74) 



The eigenstates have definite parity (even/odd) and are 
given by linear combinations of even/odd charge states. 
The CPB operating point is, without loss of generality, 
at n g — 1/2. Near this operating point, the CPB is well 
described by the reduced Hamiltonian 



H, 



CPB — 



( E c {2n g f 

E c (2n g - 1) 



-Ej/2 







-Ej/2 


E c (2n g - 2) 2 



(75) 

The reduced CPB Hamiltonian has a single odd eigen- 
state, the \q = 1) charge state, 

\o,0;n g ) = \l), (76) 

with n s -dependent eigenenergy 

E (n g ) = E c {2n g - l) 2 , (77) 

and two even eigenstates, \e,±;n g ), with energies 

E ± {n g )=Ec+E (n g )±^ 10 (n g ). (78) 

The qubit frequency depends on the gate voltage as 

wio(%) = \/^Ec) 2 (2n g - 1) 2 +E 2 j. (79) 

Note that at the operating point we have wio(l/2) = Ej 
and that the frequency rises quickly at a narrow distance 
from the optimal point, more than doubling for \n g — 
1/2 1 ~ Ej/Ec < 1. In terms of the charge states, the 
two even eigenstates are 



\e,-;n g ) = cosft|0) + sinft|2) , 
\e,+;n g ) = sin0|O) - cos6»|2) , 



(80) 



FIG. 3: Left panel: spectrum of the reduced CPB Hamil- 
tonian, Eq. (751, around the operating point n g = 1/ 2 fo r 
Ej — O.lEc- Dashed line: energy of the odd state, Eq. (77 1. 



Solid lines: energies of ground (bottom) and excited (top) 
even states, Eq. ( |78| |. Right panel: in the presence of a small 
inductive energy El, the CPB bands act as potentials in the 
quasimomentum space, see Ref. Q3] Dense horizontal lines 
represent a few energy levels near the edges of the bands. 



where 



cos ft = 



1 / 1 4E c (2n g -l) 
wio(n fl ) 



V2 



(81) 



The non- vanishing matrix elements of sin ^3/2 can be 
readily obtained using the charge basis form of this op- 
erator 



2i ^ 



|« + l)<«|-|g)<? + l 



For the states in Eqs. (76) and (80 1 we find 



[o, 0;n g \ sin^|e, ±;n g ) 



1± 



Wlo(%) 



(82) 



(83) 



We stress that the transitions are not between the qubit 
(i.e., even) states, but between the even and odd states; 
the corresponding transition frequencies are w±(n g ) — 
E c ±uj 10 (n g )/2, see Eqs. ([77] and ([78]). Therefore the 



tunneling of a quasiparticle into thcCPB changes the 
parity of the state, an effect known as "qua sipa rticle 
poisoning". 6 Substituting the matrix element (83) into 
Eq. (32) and using the high-frequency expression (41) we 
find 



r 



e,H — >o.O 



1 



Ej 
UJw(n g ) 



2E, 



2A 



(84) 



V u+{n g ) 

for the transition between even excited and odd states. 



In thermal equilibrium with T <C uj+(n g ), using Eq. (21 ) 
we obtain 



e,H — >o,0 



E, 



to w (n g ) 



AEj 



T 



uj+(n g ) 



-A/T 



(85) 
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Within our approximations, this expression reproduces 
(after implementing the corrections described in Ref. |22] 
and up to a numerical prefactor) the decay rate calcu- 
lated in Ref. [5] for the "open" qubit at the operating 
point n g = 1/2. For the transition between even ground 
and odd states the matrix element in Eq. ( 83 ) vanishes 



at the operating point. This vanishing is a consequence 
of the low-energy approximation that lead to Eq. ([6| : as 
the results of Rcfs. 5 18 show, the contributions that we 
neglect cause a finite transition rate, which is suppressed 
by a small factor of order Ec/2A in comparison with the 
transition rate from even excited to odd state. 

We note that while in all the above expressions the 
distance |2n ff — 1| from the operating point can be large 
compared to the small parameter Ej/Ec <C 1, the de- 
scription based on Eq. ( 
can be neglected, whic 
\2n 



75 1 is valid if other charge states 



g -l\< 1/2 (with \2n g 



r limits the range of validity to 
1| - 1/2 » Ej/Ec). For ex- 
ample, at 2n g — 1 ~ 1/2 the charge states |0) and |3) are 
nearly degenerate and we can expect an enhanced tran- 
sition rate r ej _| in comparison to the rate T e ^ >0 ,q 

that we have considered above. 

Finally, let us comment on the relationship between the 
transition rate in the CPB and in the inductively shunted 
Josephson junction with large charging energy [see the 
paragraph containing Eq. (59)]. As shown schematically 
in the right panel of Fig. ]3]and discussed in detail in 
Ref. HU the spectra of the two systems are distinct even 
in the limit of small inductive energy El', in the CPB 
(Ei = 0) the energy levels form bands as n g varies, 
while for any non-zero El the gate voltage n g can be 
"gauged away" and the spectrum consists of discrete lev- 
els that become denser as El decreases. Despite these 
differences, the ac responses of the two systems due to 
charge coupling agree in this limit Similarly, we now 
show agreement for the quasiparticle transition rates. We 
note that when taking the limit El — > 0, the condition 
Ej < ulc for the validity of Eq. 



for the rate r„^o 
requires that we also take Ej — > Of^TVIoreover, since the 
final state considered in deriving the rate T n ^,o is the 
lowest possible state, the corresponding final state in the 
CPB is either the even ground state at n g — or the 
odd ground state at n g — 1/2. Indeed, the width of the 
ground state (in quasimomentum space see Fig. 3 and 
Ref.[Hl) is cx (E L /E C ) 1/4 , so that as E L -> the state is 
localized at the bottom of the band. Note that following 
the same procedure detailed above it is straightforward 
to show that the transition rate r o>+ _ 



, ej o at n g = coin- 



cides with r e ^ 5. o at n g = 1/2; hence for our purposes 

the two possibilities are equivalent. At finite El, the 
total transition rate to the ground state is obtained by 



summing Eq. (591 over all initial levels n. Due to the 



Gaussian factor in the second line of Eq. (59), the num- 



ber of levels that contribute to the total rate is approxi- 
mately tJtm> cx (Ec/El ) 1 / 4 , which grows as the inductive 
energy diminishes. However, the energy of the contribut- 
ing levels tends to the charging energy, as can be seen 
by rewriting identically the argument in the exponen- 



60 
50 
40 
30 
20 
10 

-10 




FIG. 4: Potential energy (in units of El) for a flux qubit 
biased at <£> e = <E>o/2 with Ej/El = 10. The horizontal lines 
represent the two lowest energy levels, with energy difference 
e given in Eq. ( 100 1. 



tial of the Gaussian factor as — (E n — E c ) 2 /2E C ^LC: 
where E n = ululc', this agrees with frequency for the 
e, H > o, transition at n g = 1/2 in the CPB being ap- 



proximately Ec in the small Ej limit. Using Eq. (59), 



performing the sum over levels, and taking the limit 
El — > 0, we find 



lim VT n _>.o 



AEj 




(86) 



which coincides with the leading term of Eq. ( 84 ) in the 
limit Ej — > at the operating point n g — 1/2. 



B. Flux qubit 

As a second example of a strongly anharmonic system, 
we consider here a flux qubit, i.e., in Eq. ^ we assume 
Ej > El and take the external flux to be close to half 
the flux quantum, $ e « <&o/2. Then the potential has a 
double- well shape and the flux qubit ground states |— ) 
and excited state |+) are the lowest tunnel-split eigen- 
states in this potential)*! see Fig. |U The non-linear na- 
ture of the sin ip/2 qubit-quasiparticle coupling in Eq. ^ 

has a striking effect on the transition rate T_| which 

vanishes at <E> e = &o/2 due to destructive interference: 
for flux biased at half the flux quantum the qubit states 
|— ), |+) are respectively symmetric and antisymmetric 
around ip = tt, while the potential in Eq. ^ and the func- 
tion simp/ 2 in Eq. (32 ) are symmetric. Note that the lat- 



ter symmetry and its consequences are absent in the envi- 
ronmental approach in which a linear phase-quasiparticle 
coupling is assumed. 

Analytic evaluation of the matrix element determin- 
ing the transition rate [Eq. (32)] at finite $ e — 3>o/2 
is possible when Ec *C Ej and the tunnel splitting 
e is small compared to inductive and plasma energies, 
e <C 2it 2 El <C uj p ; an estimate for the splitting is given 
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below in Eq. (100). With the above assumptions we can with 



use a tight-binding approach. Neglecting tunneling the 
wavef unctions \m) Eire , £ls £l first approximation, ground 
state wavefunctions of the harmonic oscillator with fre- 
quency uip and oscillator length I = -J 8Ec /u p localized 
around the (flux-dependent) minima tp m of the potential 
energy, 



(tp\m) 



1 



1/4 



(87) 



The minima are found by solving Eq. ([49]) approximately, 
using the condition El <C Ej (which follows from the 
above assumptions) to get 



2tt 



El 

Ej 



(m-/) 



/ 



$ 



The energies of the localized states are (up to a constant 
term) 



E m = 2TT 2 E L (m - fy 



where 



E, 



E, 



El 

E L 



(89) 



(90) 



takes into account corrections small in 1/(3 <C 1. The 
above results are valid for \uiEl/Ej\ <C 1. Still ne- 
glecting tunneling, the matrix element of sintp/2 between 
states localized in different wells vanishes, but the diago- 
nal matrix clement is finite due to the shift of the minima 



away from 2-ixm, see Eq. (88 1. Using the states in Eq. ( 87 ) 
we obtain 

<j|sm||m) ~ -(-l)"v||( m - f)S mJ . (91) 

To include the effect of tunneling we allow for the pos- 
sibility of transitions between neighboring wells with am- 
plitude e/2. As we are interested in the two lowest eigen- 
states for / near 1/2, we consider only the m = 0, 1 wells 
and the effective Hamiltonian has the form 



H = 



2n 2 E L f 2 -e/2 
-e/2 2^E L (\-]f 



The eigenenergies are [cf. Eqs. |78[)-(81 )] 

E±(f) = y£ i [l + (2/-l) 2 ]±^ 10 (/) 
with the flux-dependent qubit frequency 



u w (j) = \J?+[{2nyE L {f-l/2)Y 
while the eigenstates are 



cos0|O) 
sin0|O) • 



-sin0|l) 
cos 011) 



(92) 



(93) 



(94) 



(95) 



cos 6 




(2^)2^(7-1/2) 



wio(/) 



(96) 



The tunnel splitting e entering in the above formulas can 
be estimated by noting that due to the assumption /3 3> 1 
the wells are nearly symmetric. Neglecting the asymme- 
try [i.e., considering the potential in Eq. ^ at / = 1/2], 
the width and height of the tunnel barrier are approxi- 
mately 27r(l — 1//3) and 2Ej, respectively, with 



E.j — Ei 



(97) 



To account for the height and width at El ^ 0, we treat 
the two wells as cosine potentials with renormalized coef- 
ficients. That is, we consider e ach well to be described by 
the Hamiltonian given in Eq. ( B3 ) with the substitutions 
Ej — !• Ej and Eq Ec, where 



Eq = Eq 



1 



(i - i/PY 



(98) 



Then we can use the known asymptotic formula ^ * 14 * 24 ! for 
the splitting eo in the periodic cosine potential (i.e., for 
El = 0; see Appendix [B] for a derivation of this formula) 

eo =4^(^yV^7^ (99) 



7T P \E C 



to find 



1/4 



-^d,Ej/E c 



(100) 



Here the numerical prefactor is smaller by factor of 2 in 



comparison with Eq. ( 99 ) to account for tunneling being 
;her than in a periodic potential! 2 ^ 



between two wells rat 

Turning now to the matrix element (j\ smip/2\m), the 
diagonal elements j = m = 0, 1 arc still approximately 
given by Eq. (91). Tunneling introduces finite but ex- 



ponentially small off-diagonal elements which, similarly 
to the splitting, can be calculated using the semiclassi- 
cal approximation. Using the wavefunctions derived in 
Appendix [b] we arrive at [cf. Eq. (C6)] 



(l|sin||0) 



D 



Ej_ 
Ec 



1/3 



2V2Ej 



(101) 



with D ss 1.45, see Eq. (C7 



We can now calculate the 
matrix element of sm(p/2 between qubit states |±) in 



Eq. (95) using Eqs. (91) and (101 1 to obtain 



(-|sin^!+) = tt(/-1/2) 



wio(/) 



El 
Ej 



V2n D EL 
Ej 



Ej_ 

E C 



1/3' 



(102) 
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Here the first term in square brackets is the combina- 



tion of the two intrawell contributions [Eq. (91)] while 



the second one originates from the under-barrier tun- 
neling [Eq. (101)]. Comparing Eq. (102) to numerical 



calculations, we find that near half the flux quantum, 
|/ — 1/2 1 < c/2tt 2 Ei j1 the two approaches give the same 
dependence on flux and agree on the order of magnitude 
of the matrix element, with Eq. ( 102 ) providing a smaller 



estimate than the numerics by a factor of about 2/3 . For 
|/ - 1/2| > e/(2n) 2 E L the flux dependence in Eq. pi 
via the factor (/ — l/2)/w 10 (/) can be neglected and the 
right hand side reduces to a flux-independent constant. 
However, this behavior is an artifact of our approxima- 
tions: for these larger deviations of flux from half the flux 
quantum the matrix element acquires additional flux de- 
pendence, beyond that given in Eq. p2] , once the asym- 
metry of the potential is taken into account. Moreover, 
for very small flux, |/| < (e/4:V2Tr 2 E L ) 2 , mixing of the 
sate localized in well m = 1 with that localized in well 
m = — 1 cannot be neglected and the matrix element has 
a narrow peak around zero flux. Substituting Eq. ( 102 ) 



into Eq. ( 32 ) , keeping the leading contribution, and us- 



ing the relation ( 40 1 , we find for the transition rate in 
the high-frequency regime^ 



uw l 

7T 9K 



AttEj 



1 



El 

Ec 



2/3 



(103) 



18) 



with ReY^J of Eq. 

The rate in Eq. ( lo3 1 depends on reduced flux / via 
the qubit frequency, see Eq. (94). In particular, for 
external flux equaling half the flux quantum we have 
wio(l/2) = e and the transition rate vanishes, as dis- 
cussed above. In the previous section we mentioned in 
the text after Eq. ( 85 1 that for the Cooper pair box the 



vanishing of the rate at the operating point is valid up to 
small corrections, being a consequence of the low-energy 
approximation for the tunneling Hamiltonian in Eq. (|6|. 
The same is true for the flux qubit; in the present case, 
the parameter suppressing these corrections is exponen- 
tially small, being given by e/2A. Note that if keep- 
ing in Eq. ([5| the contributions beyond the low energy 
approximation, the operators accounting for the qubit- 
quasiparticle interaction cannot be reduced to sin ip/2; 
therefore, for these additional contributions the symme- 
try argument given at the beginning of this section for 
the vanishing of the transition rate at / = 1/2 does not 
hold. 



MULTIPLE-JUNCTION QUBITS: GENERAL 
THEORY AND APPLICATIONS 



In this section we generalize the theory of Sec. |TT] to 
the case of systems containing multiple junctions. This 



generalization will enables us to consider the flux depen- 
dence of the transition rates in the two-junction split 
transmon and in the many-junction fluxonium. These 
two qubits are particular examples of the general case 
in which M + 1 junctions separate M + 1 superconduct- 
ing islands forming a loop. We use the convention that 
junction j = 0, . . . , M is between islands j and j + 1 
and identify island j = M + 1 with island j — - see 
Fig. [5] When the loop inductive energy is much larger 
than the Josephson energies of the junctions (i.e., the 
loop inductance is small), the phases are subject to the 
flux quantization constraint 



A I 



Y^fj = 27T$ e /$ . 



(104) 



3=0 



This constraint must be taken into account to derive the 
Hamiltonian H^y of the M independent phase degrees 
of freedom 6, (fc = 1, ... ,M — 1) starting from the 
Lagrangian 26 for the M + 1 constrained phases ipj 



A I 



£ M = E 



3=0 



Ejj cos tfj 



(105) 



where the dot denotes derivative with respect to time, Cj 
is the capacitance of junction j, and Ejj its Josephson 
energy. In Appendix [F] we derive the Hamiltonian as- 
suming M of the M + 1 junctions to be identical, which 
is relevant for both the split transmon (M = 1) and the 
fluxonium (M ^> 1). Explicit expressions for the Hamil- 
tonian in these two cases are presented below. 

The total Hamiltonian H of the system consist of three 
terms, as in Eq. |l]): 



H = H. 



(106) 



In addition to discussed above, the second contri- 

bution is the quasiparticle Hamiltonian 



E 

=0 



H 3 
qp ' 



qp 



EC? Q/jt aj 



(107) 



Here the index j denotes the superconducting island; 
other symbols have the same meaning as in Eq. Q and 
we assume equal gaps in all islands, A J = A. The final 
contribution to H is the tunnel Hamiltonian, given by 
the following sum [cf. Eq. Q] 



M 



H.c. 



(108) 



The transition rate between qubit states can again be 
calculated using Fermi's golden rule as in Eq. (31). We 



assume that the quasiparticle distribution functions are 
the same in all islands and that tunneling across each 
junction is not correlated with tunneling in nearby junc- 
tions - this is a good assumption if the mean level spac- 
ing in the finite size superconductors is small compared 
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E.jn Cn 



E.n Ci 




tances by 



FIG. 5: Left: schematic representation of the split transmon 
with two (possibly different) junctions. Right: in the fluxo- 
nium a weaker junction (j = 0) is connected to a large junc- 
tion array (j = 1, . . . , M). 



to the gap. Then the total rate for the transition between 
eigenstates of Hamiltonian H^y is 



M 

3=0 



(/«}l sin yl*w) 



EjjS^faf), (109) 



where for convenience we have extracted the Joseph- 
son energy prefactor from the spectral density, S qp = 
Sqp/Ej, with Sqp defined in Eq. (33 1. Similarly, the cor- 



rection 8Ei to the energy of state is given by sums 
over junctions which generalize Eqs. (45) and (46), 



5Ei = 5E iiEj + SEiw , (110) 

M 

SE i}Ej =^2Ej j {i {lj ,y\co80 j \i w )(x q p + 2x^ p ), (111) 

3=0 

M 2 

F qp (w ifc ),(112) 



3=0 



(fc w |sin-^-|i w ) 



In the next subsections we use 
Eq. (jl09h to calculate the transition rates for the split 



wher e F q p = F qp /Ej. 



transmon and the fluxonium, and Eq. (110) to find the 



frequency shift in the split transmon. The flux-dependent 
transition rate between the two lowest even and odd 
states of a split Cooper pair box has been recently con- 
sidered in Ref. [27|for gate voltage tuned at the operating 
point. 



A. Split transmon 

A split transmon consists of two junctions, j = 0, 1, in 
a superconducting loop, see Fig. [5j Therefore, there is 
only M = 1 degree of freedom, which we denote with </>, 
governed by the Hamiltonian 

i? = AE C N 2 - Ejo cos(0 - 2vr/) - E, n cos , (113) 

see Appendix[Fj Here N = -id/dcf), f = $ e /$ , and the 
charging energy Eq is related to the junctions' capaci- 



E n = 



2(C + Ci) ' 



(114) 



Note that the Hamiltonian is periodic in / with period 
1 , so we can assume | / 1 < 1/2 without loss of generality 
(i.e., we can measure the normalized flux from the nearest 
integer). After shifting <fi — > <f> + irf, the sum of the two 
Josephson terms can be rewritten as 

E J0 cos{(j)-2TTf) + Ej 1 cos<j)^ Ej(f)cos((f>-#) , (115) 

where the effective Josephson energy Ej is modulated by 
the external flux 

Ej(f) = (Ejo + E n ) cos^/)^ + d? tan 2 (^/) (116) 

with 

, Ejo — Eji 



E 



.70 



E 



.71 



and 



tani? = <itan(7r/). 
After a further shift <f> — > <fi + i} we arrive at 
= 4E C N 2 - Ej(f) cos0, 



(117) 



(118) 



(119) 



which has the same form of the Hamiltonian for the single 
junction transmon [i.e., Eq. ^ with El = 0] but with 
a flux-dependent Josephson energy, Eq. ( |116[ ). There- 
fore the spectrum follows directly from that of the single 
junction transmon (see Fig. [2]) and consists of nearly de- 
generate and well separated states. The energy difference 
between well separated states is approximately given by 
the flux-dependent frequency [cf. Eq. ( 56 )] 



w p (/) = ^%E c Ej{f). (120) 

Note that for the system to be in the transmon regime 

Ej{f) » E c (121) 

at some flux, a necessary condition is 

Ejo + En » E c . (122) 

Then we can distinguish two cases. First, in the nearly 
symmetric case of junctions with compara ble Jo sephson 
energies, \Ejq — Eji\ < Ec, the condition (121) is satis- 
fied not too close to half the flux quantum, 



\f\- 1/2 ^Ec/niEjo+Ejt) 



(123) 



On the other hand, if the Josephson energies a re su ffi- 
ciently different, \E J0 - E JX \ > E c , then Eq. ( |121| ) is 
satisfied at any flux. 
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The transition rate Ti^o between the qubit states |0), 
1 1) can be calculated using Eq. ( 109 1 if we know the rela- 
tion between tpj and <^>; the same relation is also needed 
to calculate the transition rate r o _j. e between nearly de- 
generate states - see Appendix|C 2|for details. According 



2.5x10° 



to Appendix |Fj for the variable <j> in Eq. (113) we have 
ipi = <f> and ifo — 2nf — 4>. Accounting for the two 



changes of variables performed to arrive at Eq. (119) we 
obtain 



<Po = nf - ■& 



(124) 



derivation of Eq. ( 53 1 to find 

E c 



In the transmon regime (121), we proceed as in the 

1 + cos(tt/ ± 1?) 



|sin^|i; 



w P (/) 



(125) 



where the upper (lower) sign is to be used for j = 1 
(j = 0) . Substituting this result into Eq. ( |109 ) and using 
Eq. (41) with w = 0J p (f), we find in the high frequency 
regime [cf. Eq. (55)] 



r^ = 




(126) 



For the transition quality factor we consider, as in 
Sec. |III A[ the coexistence of equilibrium and non- 
equilibrium quasiparticles [see Eq. (61)] to find 



1 

Qio 



+ 4e~ A/T cosh 




(127) 



In Fig. [6] we show with solid lines the quality factor as 
function of temperature for 4 different values of flux / 
in a symmetric transmon (d=0). As we discussed in 
Sec. Ill A an extrinsic relaxation mechanism could be 



limiting the low temperature quality factor. Character- 
izing this mechanism by a constant quality factor Q GX t 
and assuming that only equilibrium quasiparticles are 
present, the transition quality factor has the form 



1 



no, tot 




(128) 



The dashed lines in Fig. [6] show Qio,tot as a function of 
temperature for the same values of flux; the quality factor 
Qcxt is chosen so that the zero-flux curve coincides with 
the zero flux-curve described by Eq. ( |127 ). The change 
of quality factor with flux is markedly different in the 
two limiting cases (namely, presence of non-equilibrium 
quasiparticles and no extrinsic relaxation mechanism vs. 
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FIG. 6: Quality factor as function of 2T/u p (0) in a symmetric 
split transmon. Solid lines are obtained from Eq. (1271 using 
a small nonequilibrium quasiparticle density x ne = 3.8 x 1CP 7 
and a gap value such that A/cj p (0) = 6.9 (these parameters 
are taken from experiments on single junction transmorP^). 
Flux increases from top to bottom; we show curves for f — 0, 
0.2, 0.3, and 0.4, respectively. We plot Eq. dl28b with dashed 



lines for the same values of flux. The extrinsic quality factor 
is chosen so that solid and dashed lines match at / = 0. 



extrinsic relaxation with no non-equilibrium quasiparti- 
cles) described by Eqs. ( |127[ ) and ( |128[ ). Therefore the 
measurement of the temperature and flux dependencies 
of the quality factor should give indications on the pres- 
ence of non-equilibrium quasiparticles. For example, the 
low-temperature measurements reported in Ref. [25J are 
compatible with a flux-independent quality factor; to ex- 



plain the data with Eq. (127) rather than Eq. (128) one 



would need to assume a quasiparticle density that de- 
creases with increasing flux. Since magnetic fields are 
known to break pairs and thus increase the quasiparticle 
density, for the transmons considered in Ref. |28j it is un- 
likely that non-equilibrium quasiparticles are the source 
of the low-temperature qubit decay. 

The frequency shift for the split transmon is obtained, 
as in Sec. |IIIB"l by calculating the difference between cor- 
rection to energies of nearby levels, Eq. (67). The ma- 
trix elements appearing in Eqs. (Ill )-( 112 ) are given by 
Eqs. (53) and ( |69[ ) with ujxq = oj p (f), ipo = 9 + nf for 
j = 1, and ip ^d-irf for j ; = [cf. Eq. |l24|-(p5|]. 
Using those expressions we find 



SE i+ i iEj ~ 5E. hEj - (x qp + 2^4) (129) 



and 



1 ^(0)+^(/) 



16 (130) 



Then using the relation (48) and Eq. (29), we arrive at 



15 



the high frequency result 



M/) 




(131) 



At zero flux this expression agrees with Eq. (|73j) applied 
to a single junction transmon (wio = lj p , ipo — 0). How- 
ever, similarly to the flux qubit, at finite flux the split 
transmon frequency shift is sensitive to the occupation 
of the Andreev bound states, see the first term in curly 
brackets. 



B. Fluxonium 

In the fluxonium, an array of many identical junctions 
(M 3> 1) of Josephson energy Ej% ^> Eci is connected 
to a weaker junction with Ej < Ej\. Then the Hamil- 
tonian H{<f,} for the M independent degrees of freedom 
can be approximately separated into independent terms 
for the qubit phase <f> and the M — 1 phases tf>k , 



Af-i 



H 



where 



M = 
Ha = 
Hk = 

Er = 



Ha 



(132) 



fc=i 

Y2 



4E c N'-Ej cos0+ -E L ( 



±E C iNl + -Ej 1 ^ 



k > 



Eji 

M ' 



1 

Ec 



1 1 

1 . 

Eco MEci 



(133) 



see Appendix [F] There, we also give [Eq. (F5|] the rela- 
tion between the M+l (constrained) ifj variables and the 
M independent <f>k variables. Accounting for the changes 
of variables that bring the Hamiltonian in the form given 
above, we have schematically 



L j ({<p k }) + (2Trf-cf>)/M t 



,M. 
(134) 



Here Lj ({</>&}) denote linear combinations of the vari- 
ables 4>k, k = 1, . . . ,M — 1, whose specific form can be 
found in Appendix [F] but is not needed here, while we 
show explicitly the dependence of the constrained vari- 
ables ifj on the qubit phase <j>. As in the previous section, 
we take |/| < 1/2 without loss of generality. 

As an example of the calculation of the transition rate 
for such a system, we assume that the plasma frequency 
ujpi = y/8E~ciE~Ji of the array junctions is larger than 
the other relevant energy scales (namely, quasiparticle 
energy 5E and qubit frequency ujiq). Then we can take 



the many-body state of the system in the product 

form 



M-l 



(135) 



k=l 



where IV^) is a low-energy eigenstate of H<f, and |0fe) is 
the ground state wave function of the fcth osc illato r. The 
approximations used to derive in Eq. (132) imply 



that in the formula ( 109 1 for the transition rate we can 
linearize the sine for j = 1, . . . , M. Therefore, for the 



transition rate between two states of the form ( 135 ) wc 
obtain 



Fi->f - S qp (ujif) 



E 



JO 



(/</>! sin 



+E L 



(136) 



In the weak tunneling limit e <C 2-7r 2 £x -C w p = 
^/8E~cEjo (with e the tunnel splitting of the qubit states 
at / = 1/2), we can use directly the results of Sec. IV B 
the flux-dependent qubit frequency wio(/) is given by 
Eq. (|94| and the first excited state = |+) and ground 
are the linear combination of states lo- 
For the first term in 



state"17V) = 

calized in wells m = 0, 1 in Eq. 



(95) 



square brackets in Eq. (136), the matrix element is given 

i t< 



by Eq. (102). To evaluate the second term in the same 

|m), |j) - that is, states 
we have 



regime, we note that for states 
localized in wells m and j as in Eq 



(87) 



01 2 M = n 



m 1 



El_ 
E 



JO 



E L 
E 



JO 



Therefore for the states in Eq. ( 95 ) we find 

2 _ 2 



Yi-^Y — 

2 7 V EjoJ w? (/) 



(137) 



(138) 



In contrast to the matrix element of siny>/2 considered 
in Sec. |IVB| the contribution due to tunneling can be 
neglected in this case. Substituting this result and the 
leadi ng t erm from Eq. ( 102 ) into the square brackets of 
Eq. (fl36b we getP 



4 dW) 



Er 



Er 



E 



JO 



f- 



(V2d\ 



1 - 



Ejq 
E c 

E L 
Ejo 



2/3 



(139) 



In this expression, the first term in square brackets orig- 
inates from the weak junction and the second one from 
the array. Note that when considering flux near half the 
flux quantum we can neglect the first term in comparison 
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to the second and the losses due to the array dominate 
over those due to the weak junction. Keeping only the 
leading contribution in Eq. (139) and using Eq. (41) in 



Eq. (136), we arrive at the expression for the rate in the 



high-frequency regime 



/ 2A 
wio(/) 



2-kEt 



)(l/2) 
(/) 



(140) 



10 



with wio(/) defined in Eq. (94). Note that since the 



frequency increases as the reduced flux / moves away 
from 1/2, the transition rate is the largest at half the 
flux quantum. 



VI. SUMMARY 

In this work we have presented in detail a general ap- 
proach to study the effects of quasiparticles on relaxation 
and frequency of superconducting qubits. The theory 
is applicable to any qubit - the case of single-junction 
systems is considered in Sec. [TT] and the generalization 
to multi-junction ones is given in Sec. [V] Our analysis 
is valid for both thermal equilibrium quaisparticles and 
arbitrary non-equilibrium distributions, so long as the 
quasiparticle energy is small compared to the qubit fre- 
quency - this condition, not necessary in thermal equi- 
librium, ensures that quasiparticles primarily cause re- 
laxation and not excitation of the qubit. 

For single-junction qubits, we have studied in Sec. |III| 
the weakly anharmonic limit. For small phase fluctua- 
tions, both quality factor (Sec. Ill A) and frequency shift 
(Sec. Ill B ) are determined by transitions between neigh- 
boring qubit levels and can be related to real and imagi- 
nary part of the "classical" junction admittance, respec- 
tively. The small fluctuation case applies to phase and 
transmon qubits and our results in Eqs. (62) and ([73 1) 
have been successfully tested in recent experiment d 1 1 
with these qubits. For strong anharmonicity, we have 



presented in Sec. IV results for the quasiparticle transi- 
tion rate in the Cooper pair box and the flux qubit. 

We have considered two examples of multi-junction 
qubits, the two-junction split transmon in Sec. V A and 



the many-junction fluxonium in Sec. VB In particular, 
we argue that measuring the temperature and flux depen- 
dencies of the quality factor of a split transmon could help 
resolve the question of wether non-equilibrium quasipar- 
ticles are present at low temperatures, see Eqs. (127)- 
p8| and Fig. |] 
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Appendix A: Correction to energy levels 

To calculate the correction to the energy levels as pre- 
sented in Sec. |II C[ we must account for both quasiparti- 
cle and pair tunneling. Note that due to energy conser- 
vation the latter does not affect the transition rate 
between states \i) and |/) so long as Vif < 2A; for this 
reason the pair tunneling Hamiltonian H T was neglected 
in Eq. ([!]) . More generally, the total Hamiltonian of the 
single junction system is 



with 



H n + Hj 



Ho = H„ 



H T 



(Al) 



(A2) 



The Hamiltonians H v , H qp , and Hj> are defined in 
Eqs. and (|5j), respectively and the pair tunneling 

term is 



4- I e~ i iv R u L + e i ^v L u R I a R , a L ^ 



(A3) 



(L o R). 



The last term in Eq. (Al) 



H E j = Ej cos <p , 



(A4) 



is necessary to avoid "double counting" : the Josephson 
energy originates from pair tunneling.so its inclusion in 
the effective Hamiltonian H v , Eq. (|2|, must be com- 
pensated for by subtracting the same term here. We 
will show below that this treatment is justified for small 
quasiparticle density. 

In both the quasiparticle tunneling Hamiltonian Ht, 
Eq. ([5]), and the pair tunneling one in Eq. (A3), using 
the definitions given after Eq. ^ the (real) Bogoliubov 
amplitudes are 



K) 2 = i-«) 2 = ^(i + f), j 



= L,R. (A5) 



As in the main text, we assume equal gaps and dis- 
tribution functions in the leads, A L = A R = A and 
/ = f R = f. Moreover, we neglect the contributions 
of the charge mode /q(c) = (/(£) — /(— £))/2, since they 
are suppressed by the small factor SE/A <C 1 compared 
to the leading contributions due to the energy mode fs, 
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Eq. (14 1; for simplicity, in this Appendix we drop the 



subscript E. 

We want to evaluate the correction SEi to the energy 
of level i of the qubit at second order in the tunneling 
amplitude t for small quasiparticle density. Thus, Hq in 
Eq. (A2) is the unperturbed Hamiltonian, and we distin- 



we add and subtract the term obtained by setting Wik — 
in the denominator; more precisely, we define 



P 



1 



= -x lim 



1 



1 



e L - e R 2 ui-H>+ e L - to - e R e L +0J - e R 



(A12) 



guish three contributions to SEi, 

5E, = SEP + 5E[ 2) + 5E<; 3) , 



and separate in = 5E, 



8E;' d a finite term, 



(A6) 



6E 



7T 2 A 



caused respectively by Ht, Hj,, and Hej- Noting that 
the latter is already of second order in t, we treat it within 
first order perturbation theory to write 



P E 



OO pOO 

de L / de R 
A J A 



(A13) 



(fc|sin-|i) 



SE^ = Ej(i] cos 0\i) 



(A7) 



A+(e Ll e R ) 
xf(e L )(l-f(e R )) 



|cos-|i) 



1 



A_(e L ,e R 
1 

£i - e R 



The quasiparticle tunneling correction SE^ is ob- fr om a divergent one 



tained by second order perturbation theory, 



SE 



where 



y, J(MA} qP |#T|z,Mqp}| 2 



fc,{A} q 



Ex. 



qp 



E 



i,qp 



— Ek — Ei 



))qp, (A8) 



(A9) 



and the notation is the same as in Sec. [TTJ {ry} qp and 
{A} qp denote quasiparticle states, £\,qp and E rhCip their 
energies, and ((. . .)) qp averaging over {fy}q p . Performing 
the averaging, after lengthy but straightforward algebra 
we arrive at 



(A10) 



(fc|sin||z) 



A + (e L ,e R ) 



(fc|cos-|i) 



A_(e L ,e R ) 



f(e L )(l ~ f(e R )) (l-f(e L ))f(e R ) 



x 

e L — e _R. — ^ik e L — tR + ^ik 

where we introduced the functions 



A±(e L ,e R ) 



± 



I 2 

qp 



qp 



qp 



(All) 



A 2 



1 7T 2 A 



<L 



deL I de R 

A J A 

e R 



f(e L ) - f(e R ) 
£L ~ e R 



- (i| cos (p\i) 



(A14) 



A 



To obtain this expression we used the identities 



E 



|sin||i) 



and 



E 



(fe|sin 



(fcj cos — |i) 



(/s|cos-|i) 



= 1 



(A15) 



(l| cos (p\i) 



(A16) 

Equation (A13) for 5E\ can be further simplified 



using the relations 
1 



1 



A 



1 

e L — £R — ^ik 
1 



1 

1 



(A17) 



tL^ik 



£L — tR — (-L — (R 



( L — e R ~ ^ik 

t^ifc /A 
£L — £R — 



describing combinations of BCS densities of states. Both 
these functions and the lower integration limit depend on 
the self-consistent gap A qp ; however, since the integrand 



in Eq. (AlOl is at least linear in the distribution function 



/, we can neglect the gap suppression by quasiparticlcs 



[see Eq. (44)] and approximate A qp ~ A 



We note that the combinations of distribution func- 



tions in the last term of Eq. ( A10 ) restricts to low energies 



only one of the energy integrals, while the other integral 
is logarithmically divergent. To isolate this divergence, 



where the approximation is valid because the distribu- 
tion function restricts the integral over cl to low ener- 
gies above the gap. As discussed in Sec. the matrix 
elements of operators e ±l ^/ 2 describe the transfer of a 
single charge. For this reason, for a low-lying level i the 
main contribution in the sum over states k comes either 
from levels with energy difference u.^ ~ Ec *C A (when 
Eq is large compared to El, Ej), or from nearby levels 
(for small Ec)- In both cases we have ujik <§C A, since 
at large energy differences the matrix elements quickly 
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decrease; this is evident, for example, in the expressions 



for the matrix elements in Sec. Ill Then according to 
Eq. \A17h the term proportional to A_ is suppressed by 



the small parameter uiik /A in comparison to the leading 

( 1) f 

term in A+, and we can approximate SE\ as 



SE, 



(i)J 



(*|sin||i) 



WEj 
ir 2 A 

2 



A 



de L / de R 

A J A 



(A18) 



xf(e L )(l-f(e R )) 



Defining the function F qp by 
16Ej 



A P L d£L .L AV4-A' 



'A JA 

/(<*)(! -/(<*)) 



1 



1 



(A19) 



we arrive at the expression for the quasiparticle correc- 
tion to the energy SE itqp given in Eq. (46). 



(2) 

The treatment of the pair correction term SEi in 



Eq. ( A6 ) is similar to the above one for 5E { f ] . The pair 



correctio n is found by calculating the mat rix element of 
H%. [Eq. |A3|] rather than H T in Eq. |A8|; we find 



SE, 



(2) 



4E 



ir 2 A 



deL / de R 

Arm •'Ann 



(A20) 



(fc|sin-|i) 



A_(e L ,e R ) 
f(e L )f(e R ) 



cos — |i) 



(l_/( Ci ))(l_/( efl )) 



the levels which we neglectful Keeping only the second 
term in each square brackets, we write 



8EP' d + 5E?> 



SEf 



SEf, 



(A22) 



where, separating the terms independent of and propor- 
tional to the distribution function /, we have 



SEf 



AE f°° f°° 

{i\ cos <p\i)P / de L dej 

^ J Ann J A QO 



A,. 



1 



(A23) 



A 



qp 



A2 £L + £fi 

qp 



and 



SEf = ^-(i\ cos <p\i)P HdeL 

7T 2 A J A 



f(e L ) 



de R 



A 



£l + £r e L - e R 



(A24) 



In both expressions the integrations can be performed 
analytically [using in Eq. (A24) the definition ( A12)]. We 
obtain 



A 



qp {l\ COSLp\l) 



(A25) 



and 



SEf = 2Ejf(A)(i\cos0\i) 



(A26) 



Finally, using Eqs. (23 1, (44), and (A7) we arrive at 
8E? + SEf 



SE, 



(3) 



Ej (a 



■2x* p ) (z| cos^|i) 



(A27) 

which is the correction S Ej ^ , in Eq. (45). This result, 
together with Eqs. (A18|-(A19), concludes the derivation 



of the formulas presented in Sec. II C 



Note that in this expression there is a term independent 
of the distribution function, for which the approximation 
A qp ~ A is not applicable. Since + e R > 2A qp , re- 
peating the argument preceding Eq. ( A18 1 we can neglect 



Wjfe in the denominator and use identities (A15)-(A16) to 
obtain 



SE, 



.2) 



A A R 



£L + (R. 



A 2 

qp 



cos <p\i) 



A 2 

qp 



(A21) 
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*qp 



A 2 

qp 



A 2 
qp J 



Both in this expression and in Eq. (A14), the first term 



in square bracket does not depend on the level index i. 
Therefore, it leads to an unimportant common shift of all 



Appendix B: Gate-dependent energy splitting in the 



transmon 



The transmon low-energy spectrum is characterized by 
well separated [by the plasma frequency uj Pi Eq. (56)] 



and nearly degenerate levels whose energies, as shown 
in Fig. [2l vary periodically with the gate voltage n g . 
Here we aerive the asymptotic expression (valid at large 
Ej/Ec) for the energy splitting between the nearly de- 
generate levels. We consider first the two lowest energy 
states and then generalize the result to higher energies. 

Using the notation of Sec. [TTJ the transmon Hamilto- 
nian is 



(Bl) 



H v = AE C [N - n g ) -Ej{l + cos<p) 



Its eigengstates can be written exactly in terms of Math- 
ieu functions!^ However, since Ej 3> Ec a tight-binding 
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approach!^ can be used in which the two lowest (even 
and odd) eigenstates ^ e and are given by sums of 
localized wavef unctions, 



-in g 27rj 



(B2) 



*°(^;n g 



— p m 3 l f> _ 



7= ^ ~ 27r -?) e ~ 



in g 2-Kj e -i-irj 



where L ^> 1 is the number of sites, labeled with index 
j, and tp is the ground state of the Hamiltonian 



H = 4E C N 2 + V(<p) 



(B3) 



with 



F(^) = {-^ (1 + C0S ^' M <7r . (B4) 

V ' [ 0, \ip\ > 7T V ' 

This potential is such that J2j V(<P ~ 2n j) — ~Ej(l + 
cos ip). Note that the even (odd) state is a linear combi- 
nation of even (odd) charge eigenstates, as can be shown 
by considering the overlap of <3> e (°) with the ch arge eigen- 
statc e lnv ' 2 for arbitrary integer n [in Eq. (B2) what 



distinguish the odd state from the even one is the last 
exponential in the expression for 'i' , which changes the 
sign of the localized wavefunction at odd sites j] . 
The energy difference ui eo between the two states is 



(%°\H V \V°) 



(* e |i^|* e 



(B5) 



Using Eq. ( B2 ) , the contributions to uj eo due to products 



of wavefunctions if) localized at the same site cancel. The 
leading contribution to u eo originates from products of 
wavefunctions localized at nearby sites, 



w eo = eo cos (27m ff ) (B6) 



e = -4 / dtp ip(<p)ip(<p - 2tt)V(<p) . (B7) 



with 



To estimate the above integral, the behavior of the wave- 
function Tp near ip — n is needed; in this region a good ap- 
proximation is given by the semiclassical wavefunction^ 



Co 



ljj{tp) 



IvW) 



ex p [- JT rf< M<W] > 



a < ip < 7r 



Aq exp 



Ej 
2E C 



1 - 



(tp-n) 



where Cq and Aq are constants, 



Ej 
4£, 



c 



1 - 



2Ej 



COS ip , 



,<P > I" 

(B8) 
(B9) 



and a is the classical turning point defined by p(a) = 0. 
The constant Cq is determined by the normalization con- 
dition of the wavefunction, and Ao then follows from con- 
tinuity of the wavefunction. For states with large quan- 
tum number the semiclassical approximation can be used 



also in the classically accessible region \ip\ < a; the corre- 
sponding estimate for the normalization constant, which 
we indicate with , is = yJup/AEcn - see Ref . [3D 
Here we are interested in the ground state (and more 
generally in low- lying states), for which Cqq i s known 
to underestimate the normalization factorP^EH To eval- 
uate Co we note that for \tp\ <C 7r the potential V(ip) in 
Eq. ( B4 1 is well approximated by that of the harmonic 



oscillator; therefore the semiclassical wavefunction (B8) 



should match the normalized wavefunction of the har 



monic oscillator given in Eq. (87) (with tp m = 0) in the 
region a -C tp <C it . Indeed, in this region we expand the 
cosine in Eq. (B9) and rescale variables (</> = (f>y/u! p / Ej) 
to find 



d(j>p{4>) 



(BIO) 



^p 1 - 1 - In Op + \J(p 2 - l) 



1 Ej 2 1 1 , . 

2—^ ----ln(2^K 



Using this expression, and p((p) — <f\J Ej/8Eq in the 
denominator, Eq. (B8| becomes 



eV4 {8E C \ 1/8 

m ^ Co V^{-E7) ' 



— <p Ej/2u) p 



This function matches Eq. ( 87 ) by setting 

1/4 



(Bll) 



(B12) 



The last form shows that the correct normalization factor 
is larger than the usual semiclassical estimate. 

Having found the normalization constant, we now con- 
sider the wavefunction in the region near ip = tt. There 



we can further simplify Eq. (B8) as follows: we rewrite 
the integral in the exponential in the first line of Eq. ( B8 ) 



>p(<i>) . (B13) 



Then the first integral on the right hand side is 



d<t>p{4>) 



2E 



E c 



- [E(k) - (1 - k 2 )K{k)] , (B14) 



where E and K denote the complete elliptic integrals 
with modulus k, which has the value 



1 — fe = 1 



4Ej 



(B15) 



Here we are interested in the limit k — > 1, in which the 
complete elliptic integrals behave as 



^).i + V(ml-i 

K{k) ~ In 1 . 



(B16) 
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The last integral in Eq. (B13) can be approximated as also in agreement with the expression in the literature. 



d(j>v{4>) 



Ej 
2E C 



1 



4Ej 
(ir-pf 



(tt - tp) 



24^/1 - lo p /4Ej 
(|B13|)-(|B17|) into Eq. 



(B17) 



Substituting Eqs. (|B13j)-(|B17j) into Eq. |B8|), using 
p(n) ~ \J Ej/2Ec in the square root in the denominator 
of the first line, and requiring continuity of the wavefunc- 
tion, we arrive at 



A exp{ - 



Ag exp 



Ej 
■IE, 



1 - 



4E 



< 



24y/l-u, p /4Ej 

~E~7 

2E C 



4E. 



with 



tp > IT 

(B18) 



A 



(2tt) 1 /4 



8Ej 
Ec 



(B19) 

7r can be obtained by sub- 



The wavefunction near Lp 
stituting p — > —(p in Eq. (B18). We can now proceed 
with the calculation of the integral in Eq. (|B7|). Using 



Eqs. ( |B4[ ) and ( |B18 ), expanding the potential for ip < n, 
and changing the integration variable (</? — > tt — p) we 
find 



eo ~ 2EjA Q I dipip exp 
'o 



8cj p Al 



2E C 24y/l - ujp/AEj 



8E,j 
~Ec 



1/4 



-y/8E.,/E c 



(B20) 

where, going from the first to the second line, we ne- 
glect the subleading correction originating from the de- 
nominator in the argument of the exponential. The fi- 
nal exp ression for eo agrees with the known asymptotic 
formula p 1 14 1 24 1 thus validating our approach. 

The above result can be generalized to calculate the 
splitting between nearly degenerate even/odd states of 
approximate energy nui p above the gro u nd s tate by let- 
ting ui p —> uj p (2n + 1) in Eqs. (|B8|)-(|B9|) and those 



that follow [this replacement is appropriate so long as 
ujp(n + 1/2) <C 2Ej]. Matching the semiclassical wave- 
function to the excited eigenstates of the harmonic oscil- 
lator, we find that the normalization coefficient depends 



4E C 



1/4 



1/2 



i/2 



1/2 



(B21) 

Note that C n approaches Coo as n grows. Repeating the 
above calculation, we find the energy splitting 

n/2 



e„ = eo(-l) r 



SEj 
En 



(B22) 



Appendix C: Rate of parity switching induced by 
quasiparticles in the transmon 

The spectrum of the transmon, as described in Ap- 
pendix [5] comprises both well separated and nearly 
degenerate levels of opposite parity (see also Fig. |2j. 
The leading contribution to the transition rate between 
states of different parity separated in energy by (approx- 
imately) the plasma frequency is given by Eq. ( 55 ) with 



po = and is independent of n g . Here we consider the 
quasiparticle-induced transitions between the nearly de- 
generate states 1 J fe and We first consider a single- 
junction transmon to show explicitly that the rate de- 
pends on n g and is exponentially small. Next we study 
the experimentally relevant case of a split transmon; its 
rate is qualitatively different, not displaying such expo- 
nential smallness. 



1. Single-junction transmon 



According to Eq. ( 32 ) , the quasiparticle transition rate 
r o _j. e between states VP and ^> e can be written as 



'* e |sin^l*°i 



•Sqp (k-'eo) ■ 



(CI) 



This rate depends on the gate voltage n g via the states in 
the matrix clement as well as via their energy difference 
u! eo , see Eq. ( B6 ) . For the matrix element we use Eq. ( B2 ) 
to find 



:* e |sin^|*°l 



|sin(27m 9 )| s , 



where 



dip ij){p)^i{p — 2tt) sin 



(C2) 



(C3) 



The matrix element in Eq. (C2| vanishes at half integer 



values of n g , as in the case of the Cooper pair box [see 



Eq. (83)]. In fact, the vanishing holds at arbitrary ratio 
Ej/Ec, as can be shown using the symmetry proper- 
ties of Mathieu functions. For example, at n g — 0, 1/2 
the t wo lowest eigenstates of the transmon Hamiltonian, 
Eq. (Bl I, can be written in the charge basis a*21 



* e > = E A \ 2m ) + 1 - 2m > 



m=0 
oo 

= E A ^ + i [|2m + 1} + | - (2m + 1)) 

m=0 



(C4) 
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and 



* C > = E4 1 i + i[|2m + 2) + |-2m)' 



m=0 

oo 



(C5) 



m=0 



*°)=E4! |2m + l) + |-2m + l) 



respectively, where the coefficients A^, A^l+i depend 
on the ratio Ej/Ec- 35 Using the charge basis representa- 
tion of sin p/2 in Eq. ( 82 ) it is easy to check the vanishing 
of its matrix element between the above states for both 
values of n g . 

In the transmon limit Ej/Ec 3> 1 under consid- 
eration, the product of wavefunctions localized at the 
same site does not contribute to the matrix element in 



Eq. ( C2 ) : the intrawell integral vanishes because ip 2 (tp) is 
a symmetric function (ip being the ground state of a sym- 
metric potential) which is multiplied by the antisymmet- 
ric function siny>/2; the vanishing of the intrawell term 
has thus he same origin of the vanishing of the matrix el- 
ement for a weakly anharmonic qubit at zero phase bias, 
see Eq. (53 1 with n — m and po = 0. To estimate the in- 



terwell contribution s in Eq. (|C3j), we use Eq. (B18) and 
that near ip = it we have sin Tpjl ~ 1. After changing 
integration variable (p> — > 7r — p) we arrive at 



s ~ AAq I dp exp 



E, 



2E C 24^/1- uj p /4Ej 



D 



(C6) 



where (with T denoting here the gamma function) 
D = 2 1 / 6 3~ 2/3 r i 1.45. 



(C7) 



Due to the factor eo in Eq. (C6), the transition rate in 



Eq. (CI I is indeed exponentially small. Turning now to 



the factor S qp in Eq. (CI), we note that its argument, 



w eo , is usually small due to its exponential suppression 
at large Ej/E c , see Eq. ( |B20j ). Therefore the "high fre- 
quency" condition uj eo 3> SE (with SE the characteristic 
quasiparticle energy) is in general not satisfied and use 
of Eq. (41) expressing S qp in terms of the quasiparticle 



density is not appropriate. In thermal equilibrium, one 
can use Eq. (35) for arbitrary ratio uj eo /T. Assuming 



we rewrite Eq. (CI) as 



e -C T, using Eq. (17 1, Eq. (35), and the above results, 



16 E 



J e -A/T 



In 



4T 



eo cos(27rn fl 



IE 



Ec 
Ej 



1/3 



D 



eo 



sin 2 (2Tm g ) . 



(C8) 



Generalization of this result to the transition rate r^,_> e 
between nearly degenerate states of higher energy is ob- 



tained by the substitution eo ~ ► £n- Except at the de- 
generacy points n g = 1/4, 3/4 (where this expression di- 
verges), we can estimate the rate in order of magnitude 
by assuming sin(27rn s ), cos(27m ff ) s» 1. For low-lying 

(n) 

states, this estimate shows that the rate To-te is small 
compared to the rate determining the relaxation 

time of the transmon [see Eq. (55 )]. This smallness is due 
to the exponentially suppressed o — > e matrix element, 
Eq. ( C6 1 , as function of the ratio Ej /Eq , in comparison 
with the weak power-law suppression of the 1 — > matrix 



element as given by Eq. (53) with ipo = 0, m = 1, and 



n = 0. The relationship between the two rates is qualita- 
tively different in the split transmon, as we discuss next. 



2. Split transmon 

The above calculation of the even/odd transition rate 
in the single junction transmon can be easily modified 
to yield the rate for a split transmon. As discussed in 
Sec. IV Al the effective Hamiltonian and therefore the 



form of the eigenstates are the same in the single and 
split transmon. The difference between the two cases 
arises in the evaluation of the matrix elements pertain- 
ing to each junction [cf. Eq. (125)]. For the even/odd 
matrix element we find 



I* e |sin^|*°; 



1 - cos(tt/ ± i?) 



(C9) 



where the upper (lower) sign applies to junction j = 1 
(j = 0), / is defined in Eq. @ and i9 in Eq. pl8| >. In 
contrast with the single-junction transmon case consid- 
ered above, here the matrix element is dominated by the 
intrawell contribution having the same form of Eq. ( 53 ) 
at n = m = and finite phase bias nf ± <&. Substituting 



Eq. ( C9 ) into Eq. ( 109 1 and assuming thermal equilib- 
rium quasiparticles [cf. Eq. (35)] we obtain 



8(Ej + Ejx) 



e -A/T e ^/2T Ko 



1 - 



2T 



4W\ (C10) 
"2(°) I 



with ujp(f) given in Eq. (120) and w eo in Eq. (B6). As 

(n) 

before, the rate T -ie of transitions between nearly de- 
generate levels of higher energy is obtained upon the sub- 
stitution oj Pn — > e n sin(27m g ) in Eq. (C10). Note that the 



rate vanishes at integer multiples of the flux quantum; at 
those values of flux, exponentially small contributions to 
the matrix element analogous to those calculated above 
should be included. At non-integer values of reduced flux 
/, Eq. ( C10 ) should be compared with the transition rate 
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between qubit states induced by thermal quasiparticles, 



= 8(Ejq + Eji) c - A/T( j*.(f)/xr Ku ( Mf)\ 

7T \ 2T 



M/) I ^(0) 



(Cll) 



obtained using Eq. (125). The ratio between these two 
quantities, 



2T 



£O P (/)«p(0)-w2(/) 



r 



1^0 



.(/)/2T# (feffl) E C c^(0)+w|(/) 



(C12) 

depends on temperature through the first factor on the 
right hand side. Experimentally, measurements for the 
rate are performed near n g = 1/2, so that the rele- 
vant even/odd frequencies are u eo ~ eq, e\] they are 
generally 2-3 orders of magnitude smaller than u> p (f) 
(~ 27r x 4 GHz), while the latter is usually larger than 
twice the temperature (T ~ 20 — 200 mK) . Under these 
conditions, the first factor in Eq. (C12) can be approxi- 



mated, in order of magnitude, by 5 to 10. The last factor 
in Eq. (C12 1 varies between at / = and 1 at / = 1/2; 
as flux is used to suppress the qubit frequency from its 
maximum value (> 10 GHz), we can approximate the last 
factor by 1/2. Finally, the central factor can be rewrit- 
ten as ^8Ej(f)/Ec] since Ej(f)/Ec usually is varied 
between 10 and 30, we arrive at the order-of-magnitude 
estimate 



l->0 



20-80 



(C13) 



in the experimentally relevant ranges of parameters. This 
is an example of the more general statement that, except 
close to integer values of /, the even/odd transition rate 
in a split transmon is faster than its decay rate. This 
result is qualitatively in agreement with experimental 
boun ds for the even/odd transition rate in split trans- 
monS ™ 



Appendix D: Matrix elements for the harmonic 
oscillator 



In this Appendix we present analytic expression for the 
matrix elements of sin ip/ 2 between harmonic oscillator 
states |n) and \m). Let us introduce the displacement 
operator 



a C« — M a 



(Dl) 



where a (a)) is the harmonic oscillator annihilation (cre- 
ation) operator. The matrix elements of D are^ 



(m\D(fi)\n) = < 



/' 72,/"! ( M )m- , ;/ 



n—m j ( n m ) 



n t (m-n) 



(H 2 )> 

m < n 

(ImI 2 ), 

m > n 
(D2) 

where are the generalized Laguerre polynomials. 

Since the position operator is (p — £(a + o')/v2, where 
I = 1/y/mu! is the oscillator length for an oscillator of 
mass m and frequency u), we can write 



( ii 



\2v / 2 



(D3) 



Note that for the harmonic oscillator described by 



Eq. ( 50 ) we have 




(D4) 



To allow for fluctuations around a finite phase, we shift 
tp — > tpa + <p in the argument of sine and rewrite the 
resulting expression in terms of exponentials, 



sin ■ 



fa + <P 



1 

2i 



e iip /2 e iip/2 _ e -i<p /2 e -itp/2 



(D5) 



Then using Eqs. (D2|-(D3) we find for m < n 



2 



;///! sin— -|n)=e 1 /16 



xL 



(n — m) 



Sill 




(D6) 



The matrix element for m > n is obtained by exchanging 
n f-> m in the right hand side. Equation (53) can be 



obtained from Eq. ( D6 1 by Taylor expansion for small £, 



which for the Laguerre polynomials gives 

4«\ x )= { -^^- ( {m Z a)l 1 ^ + ^ 2 )- (D7) 

m v ' m\a\ (m- l)!(a + 1)! v ' v ' 



Equation ( 58 1 follows from Eq. ( D6 ) with m = using 



Using Eq. (D2) we can also find the expectation value 



of the operator cos^. After shifting the phase variable 
as done above and since the expectation value of sine 
vanishes by symmetry, we find 



(n\ cos (ipo + tp) \n) = cos tfo( n \ cos (p\n) 



(D8) 



Writing the cosine in exponential form, using e lip = 
we arrive at 



(n| cos (tp + tp) \n) = cos tp e l2/4 L { ° ] (^j 



(D9) 
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Appendix E: Matrix elements for the transmon 



Consider now the case m = n — 1 in Eq. ( E6 ) . Using 
Eqs. pi|, |E5|), and the leading term in Eq. fE8[), the 



Here we want to show that corrections to Eq. ((53| for central term in the right hand side is approximately 



the transmon (ipo = 0) are of cubic order in EcJuJp, as 
claimed in the text following that equation. The trans- 
mon Hamiltonian is given by Eq. (Bl| and we neglect 
exponentially small corrections by setting n g — (see 
Ref. [Hand Appendices |B| and [C]) . Numbering the eigen- 
states \ip n ) starting with n — for the ground state, even 
(odd) numbered states are even (odd) functions of if, due 
to the symmetry of the potential energy. Since siny>/2 
is an odd function, the matrix element between states of 
the same parity vanishes, 

(Vn±2,-|sin||^> =0, j = 0,1,2, ... (El) 

Due to the smallness of the charging energy, Ec -C Ej , 
as a first approximation we can expand the Josephson 
energy in Eq. ( Bl I up to the fourth order in ip. In terms of 
creation/annihilation operators (cf. Appendix [P] - note 
that in the present case I = 2\[2^/ 'Ec /uj p <C 1), the 
approximate transmon Hamiltonian is 



H = H Q + dH, 
Ho = uj v [ a* a + 



5H = - E ^(a- 
12 V 



(E2) 
(E3) 

(E4) 



To first order in Ec /w p , expressed in terms of harmonic 
oscillator states the transmon eigenstates are therefore 



|^n) = |«) + \Hn) 
(A 

En — E,„ 



(E5) 



l^) = -Eli) 01 ^ |n 



E n = oj p I n 



and including the first anharmonic corrections to the 
eigenstates the matrix elements are 



(■0ml sin-IVVi 



(m\sm~\n) 



h (m| sin||(5-0 n 
-(<^ TO |sm~|n) 



(E6) 



Using Eq. ( D6 ) , we find that the leading contribution to 



the first term on the right hand side is 



(n±(2j + l)|sin§|n> ex 

2 \uj p 



i+1/2 



j = 0,1,2,.. 



(E7) 

Since we are interested in calculating the square of the 
matrix elements up to second order in Ec/u p , we can 
neglect transitions with j > 1. For j — 0, using Eq. (D6) 
at next to leading order we find 



(n ± 1| sin —\n 




(E8) 



(n - 1\ sin r |^ n ) ~ _( n _ i| sin ii\ n _ 2 ) K — 



\_ ( E, 
24 



(j^j ' V^l(n~2\(a + a^) i \n) . 

(E9) 



To calculate the last factor we note that 

(a + a)) 2 \n) = y 1 n{n - l)\n - 2) + (2n + l)\n) 
+ ^(n + l)(n + 2)\n + 2). 



(E10) 



Shifting n — > n — 2 and taking the scalar product we 
arrive at 



(n - 2| (a + a f ) 4 |n) = Vn(n-l) (n - ^) , 



(Ell) 



and substituting this expression into Eq. ( E9 1 we obtain 



(n - 1| sin^l^n 



1 (E c \ 



3/2 



6 \ uip J 
Proceeding as above we also find 



ln(n — 1) I n 



1 



2 

(E12) 



(<fVn-i|sin^|n) = \ ( — 
2 6 \ oj. 



c 



3/2 



Vn(n + 1) [ n + - 



Finally, substitution of Eqs. (E8), (E12|, and (E13I into 



Eq. (E6) gives 



(E13) 



(0„_ 1 |sin||V„) = Jn^+o(^ 

* V w p v w p 



5/2 



(E14) 



Repeating the above calculations for the case m = n+1 
and using Eq. (E7) we conclude that the square of the 
matrix element is 



(i/j m \ sin - \ip n ) 



Er 



[n8 m ,n-i + {n + l)S m , n 



+U 



-O 



Er 



Co',, 



(E15) 



Appendix F: Multi-junction Hamiltonian 

The aim of this Appendix is to derive the Hamiltonian 
for a multi-junction system starting from the Lagrangian, 
Eq. ( 105 ) . We consider a loop of M + 1 junctions and as- 
sume M of them, denoted by index j with j = 1, . . . , M , 
to be identical, so that their capacitances and Josephson 
energies are, respectively, Cj = C\ and Ejj — Eji for 
1 < jf < M. These M junctions will be referred to as the 
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array junctions to distinguish them from the j = junc- 
tion, whose capacitance Co and Josephson energy Ejq 
can differ from those of the array junctions. 

While the system comprises M+l junctions, there are 
only M independent degrees of freedom, due to the flux 
quantization constraint, Eq. (104 1. Using that equation 



to eliminate the phase ipo, the Lagrangian is 

2 



5W|> ? (F1) 



M \ M 

+E,jo cos \^2<Pj- 2vr$ e /$ + E JX ^ cos ipj . 



We introduce a new set {(/>} of M independent variables 

(F2) 



M 

E^ 



M-l 
i=i 



, k = 1,...,M- 1, (F3) 



where 



a= 1 



1 



1 



IM) M - 1 
The inverse transformation is given by 

M-l 



(F4) 



<Pk = <Pk- a 



E 

M-l 



1 



k = 1, 



,M-1, 



;=i 



1 



(F5) 



In terms of the M variables (j), (f>k (k = 1, . . . , M — 1) the 
Lagrangian is 

M-l 



l 

8e2 



fc=i 



with potential energy 
U ({<(>}) = - Ejo cos (0 - 27T$ e /$ ) 



(F6) 



M-l 



M-l 



(F7) 



- Ejx ^2 cos [^k - a ^2 c/)i + — 

k=l \ 1 = 1 , 

Introducing the M conjugate variables N — dC^/dcf) 
and Nk — dC^/dtfik (k = 1, . . . ,M — 1), the Hamiltonian 
is 

M-l 

H {<p} =N<P+J2 N k4>k - C {4>} 

k=1 (F8) 

M-l v ' 

- AE C N 2 + AE C1 N 2 k + U({</>}), 
k=l 



where 



En, = 



2(C + C 1 /M) ' 



Eci — 



2Ci 



(F9) 



The Hamiltonian in Eq. (|F8j) governs the dynamics of 
the M independent degrees of freedom of the M+l junc- 
tion system with flux quantization and M identical array 
junctions. For a two junction system we have M = 1 
and all the sums in Eqs. (F7)-(F8) are absent. Then the 



Hamiltonian is that given in Eq. (113) 



1. Fluxonium 

The fluxonium consist of M + 1 junctions such that a 
"weak" junction j = with Ejq < Ejx is connected to 
a large array of M junctions (M 3> 1) with small phase 
fluctuations, Eci <C Ej\. These conditions enable us to 
drastically simplify the last two terms of the potential 
energy U({4>}) for the M independent variables (j), 4>k 
(k = 1,...,M- 1) in Eq. pTfr . 

We consider small fluctuations of variables 4>k around 
the configuration (pf. — 0, k — 1, . . . , M — 1, which is an 
extremum of f7 for any value of [as can be checked by 
differentiating U with respect to (pk and using Eq. (F4)]. 



We further assume that typical values of <f> are small com- 
pared to 2ttM (note that since M is large, this weak re- 
striction on <p and its fluctuations still allows for phase 
slips through the weak junction). Then we can expand 
the last two terms in Eq. (F7) to quadratic order in 
and cf>/M to find 



U{{4>}) ~ - Ejo cos (0 - 27r$ e /$ ) + l -E L 



1 



-E 



M-l 

fe=i 



with 



(F10) 



(Fll) 



Hence in this approximation the Hamiltonian ( F8 ) for 
the M + l junction fluxonium separates into independent 
Hamiltonians for each of the M unconstrained variables 



M-l 



(F12) 



fc=i 



H+ = 4E C N 2 - E JQ cos(</> - 2^$ e /$ ) + ~E L <, 



H k = AE cl Nl + -Ej^l 



Up to a change of variable 



27r$ e /<I>o - and re- 



definitions of symbols, Ha, coincides with H v of Eq. &2 



The relations in Eq. ( 133 ) between the parameters of the 
M + l junctions and the energies Eq and El entering the 
effective qubit Hamiltonian follow from Eqs. (F9 1 and 
(Fill, respectively. 
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